Put-Call-Parity is the arbitrage relationship between (European) put, call, and forward.
$$PV(Call) - PV(Put) = S_0 - PV(K)$$With cash flows at maturity $T$:
$CF(T, Call) = max(S_T - K, 0) = (S_T - K)^+$
$CF(T, Put) = max(K - S_T, 0) = (K - S_T)^+$
$CF(T, Fwd) = S_T - K$
For example, if a call option is overpriced w.r.t. its replicating portfolio, i.e. Fwd + Put, sell the call and buy the replication. The resulting shift in supply and demand must adjust prices.
"""Spot price vs cash flow (put-call parity)
"""
import numpy as np
import matplotlib.pyplot as plt
# see workbook task 4 ff.
S_T = np.linspace(0,100)
K = 50
# portfolio
CF_T = np.maximum(S_T - K, 0)
CF2_T = - np.maximum(K - S_T, 0)
# the rest is plotting
plt.figure('PCP', figsize=(9,4))
plt.plot(S_T, CF_T)
plt.plot(S_T, CF2_T)
plt.axhline(y=0, linewidth=1, color="black")
plt.xlabel("Spot price")
plt.ylabel("Cash flow")
plt.annotate("K", (K,-5))
plt.show()
Heidorn, T.; Mokinski, F; Rühl, C.; Schmaltz, C.(2015): The impact of fundamental and financial traders on the term structure of oil Energy Economics 46 (2015) pp. 276 – 287
Interest rates are quoted per annum (p.a.) but speak different “languages”:
Money Market (MM):
$C_0(1+r_{mm}T_{mm})=C_T$
ISMA (International Securities Market Association - Europe)
$C_0(1+r_{isma})^T=C_T$
SIA (Securities Industry Association - US)
$C_0(1+r_{sia}/2)^{2T}=C_T$
Continous compounding (academia, financial engineering)
$C_0 e^{rT} = C_T$
T is calculated by dividing the number of days in perid T by a year basis, e.g.
$T = act/360$ (Money Markets - DE/US/...)
$T = act/act$ (Bond - ISMA)
$T = 30/360$ (Swap)
These are some examples for day count conventions, which might differ from country to country. Observe that for $act/act$ and $30/360$, a full year must always yield $T=1$.
"""Interest rate conversion (translation)
"""
import numpy as np
r_isma = 0.1
r_mm = (r_isma*360/365)
r_sia = ((1+r_isma)**0.5-1)*2
r_cont = np.log(1+r_isma)
print(f'r_isma = {r_isma:.4f}')
print(f'r_mm = {r_mm:.4f}')
print(f'r_sia = {r_sia:.4f}')
print(f'r_cont = {r_cont:.4f}')
r_isma = 0.1000 r_mm = 0.0986 r_sia = 0.0976 r_cont = 0.0953
Semi-annual compounding has value only for didactical puposes, this is, when increasing the number of compounding periods $m$
$$C_0 \left( 1+\frac{r}{m} \right)^{mT}=C_T$$"""Continous compounding
"""
import numpy as np
import matplotlib.pyplot as plt
def FV(m):
PV = 1
r = 1
T = 1
return PV * (1 + r / m)**(T * m)
#FV_m = list(map(lambda m: FV(m), range(1,250)))
FV_m = FV(np.arange(1,250))
# plotting
plt.figure('exp', figsize=(9,4))
plt.plot(FV_m)
plt.axhline(y=np.exp(1), linewidth=1, color="red")
plt.xlabel("m")
plt.ylabel("FV")
plt.show()
print(f'FV(m = 100000) = {FV(100000):.4f}')
print(f'EXP(1) = {np.exp(1):.4f}')
FV(m = 100000) = 2.7183 EXP(1) = 2.7183
To replicate a forward consider a portfolio that consists of underlying $S$ and cash $K$, which today is worth
$$f_0 = S_0 - K e^{-rT}$$and at maturity $T$
$$f_T = S_T - K$$The forward rate $F$ at initiation is set such that $f_0 = 0$, hence,
$$S_0 - F e^{-rt} = 0 \Leftrightarrow F = S e^{rT}$$For many assets, additional factors affect the forward price, such as:
Considering all factors, we have
$$f_0 = S_0 e^{(-q-y+u)T} - K e^{-rT}$$and
$$F = S e^{(r-q+u-y)T}$$"""Forward valuation
"""
def forward_price(spot, maturity, rate, income=0, storage=0, convenience=0):
return spot * np.exp((rate - income + storage - convenience) * maturity)
def forward_mtm(strike, spot, maturity, rate, income=0, storage=0, convenience=0):
return spot * np.exp((- income + storage - convenience) * maturity) - strike * np.exp(- rate * maturity)
# 12
print(f'E12: forward_mtm(55, 60, 1/2, 0.04) = {forward_mtm(55, 60, 1/2, 0.04):.2f}')
# 14
print(f'E14: forward_mtm(27, 25, 6/12, 0.1, 0.04) = {forward_mtm(27, 25, 6/12, 0.1, 0.04):.2f}')
print(f'E14: forward_price(25, 6/12, 0.1, 0.04) = {forward_price(25, 6/12, 0.1, 0.04):.2f}')
E12: forward_mtm(55, 60, 1/2, 0.04) = 6.09 E14: forward_mtm(27, 25, 6/12, 0.1, 0.04) = -1.18 E14: forward_price(25, 6/12, 0.1, 0.04) = 25.76
The risk of imperfect hedging.
Consider a long position in a risky asset $S$ to be hedged by selling (short) future $F$. We will hedge the asset in $t=1$, close the hedge in $t=2$ and analyze the basis $b$:
$b_1 := S_1 - F_1$
$b_2 := S_2 - F_2$
Value of the portfolio in $t=2$
$S_2 + F_1 - F_2 = \underbrace{F_1 + b_2}_{b_2 = S_2 - F_2}$
$b_2$ is unknown in $t=1$, hence, risky. Now also consider the underlying of the future $S^*$
$S_2 + F_1 - F_2 = F_1 + \underbrace{(S^*_2 - F_2)}_\text{time basis} + \underbrace{(S_2 - S^*_2)}_\text{underlying basis} \hspace{1cm} | S^*_2 - S^*_2 $
Reasons for basis risk:
"""Optimal hedge ratio vs linear regression beta
"""
import numpy as np
from numpy.random import randn
from numpy.polynomial.polynomial import polyfit
from scipy.linalg import cholesky
import matplotlib.pyplot as plt
corr = np.array([[1,.75],[.75,1]])
F, S = np.dot(randn(1000,2), cholesky(corr)).T
F = F * 2 # Assume future is about twice as volatile as the asset
constant, beta = polyfit(F, S, 1)
plt.figure('mvhegde')
plt.xlabel('Future (F)')
plt.ylabel('Asset (S)')
plt.scatter(F, S)
plt.plot(F, constant + beta * F, 'r-')
plt.show()
Hedge ratio for hedging asset S with future F
$$Var = \sigma^2_S + h^2\sigma^2_F-2 h \rho \sigma_S \sigma_F$$$$\frac{\delta Var}{\delta h} = 2 h \sigma^2_F - 2 \rho \sigma_S \sigma_F = 0$$$$...$$$$h = \rho \frac{\sigma_S}{\sigma_F}$$rho = np.corrcoef(F,S)[0][1]
sd_F = np.std(F)
sd_S = np.std(S)
print(f"rho(F,S) = {rho:.2f}")
print(f"Linear regression beta: {beta:.4f}")
print(f"Optimal hedge ratio h : {(rho*sd_S/sd_F):.4f}")
rho(F,S) = 0.74 Linear regression beta: 0.3800 Optimal hedge ratio h : 0.3800
A discount factor is today's price for 1 monetary unit in t, hence, a zero bond with face value 1 and maturity t.
$$C_0 = C_t DF_t \implies DF_t = \frac{C_0}{C_t}$$The discount factor does not speak a language! However, it can be translated into any interest rate language, e.g.
$$C_0 = C_t \frac{1}{(1 + r)^t} \implies DF_t = \frac{1}{(1+r)^t} \implies r = \left(\frac{1}{DF_t}\right)^{1/t} - 1$$And as usual
$$PV = \sum^{T}_{t=1}\frac{C_t}{(1+r_t)^t} = \sum^{T}_{t=1}C_t DF_t$$"""Discounting for plain vanilla coupon bonds
"""
discount_factors = [.9091, .8116, .7118] # Slide 55
def fair_price(coupon):
# sum(discount_factors) is an annuity factor
# the last discount factor [-1] discounts the face value
p = coupon * sum(discount_factors) + discount_factors[-1]
return p
print(f'Fair price 6% 3Y bond: {fair_price(0.06)*100:.4f}')
print(f'Fair price 12% 3Y bond: {fair_price(0.12)*100:.4f}')
Fair price 6% 3Y bond: 85.7750 Fair price 12% 3Y bond: 100.3700
"""Coupon (nominal yield) vs yield to maturity
"""
from scipy.optimize import minimize
import matplotlib.pyplot as plt
# see workbook task 17 ff.
#discount_factors = [.9524, .89, .8163] # Task 17
discount_factors = [.9091, .8116, .7118] # Slide 55 (normal term structure)
#discount_factors = [1/(1.12), 1/(1.10)**2, 1/(1.08)**3] # inverted term structure
zero_rates = []
for T, DF_T in enumerate(discount_factors):
T = T + 1
r_T = (1 / DF_T) ** (1 / T) - 1
zero_rates.append(r_T)
print(zero_rates)
# face value = 1!
def IRR_coupon(coupon, discount_factors):
# calculate PV from market prices
def fair_price(coupon):
p = coupon * sum(discount_factors) + discount_factors[-1]
return p
# calculate PV from irr
def IRR_price(irr, coupon):
T = len(discount_factors)
annuity_factor = (1 - (1 + irr) ** -T) / irr
p = coupon * annuity_factor + 1 / (1 + irr) ** T
return p
# target function to be minimized
def NPV2(irr):
return (fair_price(coupon) - IRR_price(irr, coupon)) ** 2
# set NPV = 0 and return both PVs and the irr
opt = minimize(NPV2, x0=0.01, bounds = ((None, None),))
return [fair_price(coupon), IRR_price(opt.x[0], coupon), opt.x[0]]
#print(IRR(.1, discount_factors))
# lets plot different coupon rates against irrs
x_coupons = range(100)
y_irrs = list(map(lambda c: IRR_coupon(c/100, discount_factors)[2]*100, x_coupons))
plt.plot(x_coupons, y_irrs)
plt.axhline(y=zero_rates[-1]*100, linestyle='--', linewidth=1, color='grey')
plt.xlabel("Coupon rate in %")
plt.ylabel("YTM (IRR) in %")
plt.show()
[0.09998900010999878, 0.11001534051800954, 0.11998964004232815]
def YTM(price, coupon, maturity):
# calculate PV from irr
def YTM_price(irr, coupon, maturity):
T = maturity
annuity_factor = (1 - (1 + irr) ** -T) / irr
p = coupon * annuity_factor + 1 / (1 + irr) ** T
return p
# target function to be minimized
def NPV2(irr):
return (price - YTM_price(irr, coupon, maturity)) ** 2
# set NPV = 0 and return both PVs and the irr
opt = minimize(NPV2, x0=0.01, bounds = ((None, None),))
return opt.x[0]
ytm_21 = YTM(100.5/100, 2/100, 5)
print(f'YTM = {ytm_21 * 100:.2f}%')
print(f'Liability swap level = Euribor {(ytm_21 - 0.02) * 10000:.0f} basis points')
ytm_22 = YTM(101.0/100, 3/100, 5)
print(f'YTM = {ytm_22 * 100:.2f}%')
print(f'Asset swap level = Euribor {(ytm_22 - 0.0205) * 10000:.0f} basis points')
YTM = 1.89% Liability swap level = Euribor -11 basis points YTM = 2.78% Asset swap level = Euribor 73 basis points
Forward rate fixed in $s$ payed in $t$ has arbitrage condition
$$(1+r_s)^s(1+r_{st})^{t-s} = (1+r_t)^t$$$$r_{st} = \left(\frac{(1+r_t)^t}{(1+r_s)^s}\right)^{1/(t-s)}-1$$Again, a discount factor $DF_t$ is today's value of 1 that is paid in $t$. Conversely, $1 / DF_t$ is a future value. This allows us to calculate "forward discount factors" by discounting $DF_t$ and compounding a shorter period $s$ with $1/DF_s$. Discounting with a forward discount factor or a forward rate must yield the same, hence,
$$\begin{aligned} \frac{DF_t}{DF_s} &= \frac{1}{(1 + r_{st})^{t-s}}\\ &... \\ r_{st} &= \left(\frac{DF_s}{DF_t}\right)^{1/(t-s)}-1\\ &= \left(\frac{(1+r_t)^t}{(1+r_s)^s}\right)^{1/(t-s)}-1 \end{aligned}$$Agreement to exchange a fixed ($C$) against a variable ($L_t$) payment, assuming annual payments
$$\begin{aligned} PV(Fix) &= C \sum_{t=1}^T DF_t\\ PV(Var) &= \sum_{t=1}^T L_t DF_t\\ &= \sum_{t=1}^T \left(\frac{DF_{t-1}}{DF_t} -1 \right) DF_t\\ &= 1 - DF_T \end{aligned}$$Price the fair swap rate $C$ such that,
$$\begin{aligned} PV(Fix) &= PV(Var)\\ C \sum_{t=1}^T DF_t &= 1 - DF_T\\ C &= \frac{1 - DF_T}{\sum_{t=1}^T DF_t} \end{aligned}$$Mark to market an existing swap with $C^*$ (the old swap rate) by comparing $PV(Fix)$ and $PV(Var)$
$$\begin{aligned} PV(Receiver) &= PV(Fix) - PV(Var)\\ &= C^* \sum_{t=1}^T DF_t - (1 - DF_T) \\ &= (C^* - C) \sum_{t=1}^T DF_t \end{aligned}$$Swap rate $C$ is also the par coupon bond rate that fulfills
$$\begin{aligned} 1 &= \sum^{T}_{t=1}C_t DF_t\\ &= C \sum^{T}_{t=1}DF_t + DF_T \implies C = \frac{1 - DF_T}{\sum^{T}_{t=1}DF_t} \end{aligned}$$Use one curve to
Bootstrapping, i.e. creating zero bonds from a market of coupon bearing instruments
$$DF_t = \frac{1 - C_t \sum_{i=1}^{t-1} DF_i}{1 + C_t}$$"""Bootstrapping from swap rates (assuming annual payments)
"""
def swap_bootstrap(swap_rates):
discount_factors = []
for s in swap_rates:
# See Slide 61
DF_T = (1 - s * sum(discount_factors)) / (1 + s)
discount_factors.append(DF_T)
zero_rates = []
for T, DF_T in enumerate(discount_factors):
# Python starts counting at 0, hence, +1
T = T + 1
# See Slide 48
r_T = (1 / DF_T)**(1 / T) - 1
zero_rates.append(round(r_T,4))
return [discount_factors, zero_rates]
example = swap_bootstrap([0.04, 0.05, 0.06])
print(example)
[[0.9615384615384615, 0.9065934065934065, 0.8376529131246111], [0.04, 0.0503, 0.0608]]
"""Single curve swap pricing
"""
discount_factors = swap_bootstrap([0.04, 0.05, 0.06])[0]
def swap_pricing(discount_factors):
# See Slide 66
s_T = (1 - discount_factors[-1]) / sum(discount_factors)
return s_T
print(f"3Y swap rate = {swap_pricing(discount_factors)*100:.2f}%")
3Y swap rate = 6.00%
"""Single curve forward swap pricing
"""
discount_factors = swap_bootstrap([0.04, 0.05, 0.06])[0]
# Please note that discount_factors is a vector that starts with the first discount factor (t = 1)
s_13 = (discount_factors[0] - discount_factors[-1]) / sum(discount_factors[1::])
print(f"1X2 swap rate = {s_13*100:.2f}%")
1X2 swap rate = 7.10%
Since the financial crisis exists a considerable spread between 6M Euribor and 6M Eonia Swap
Since there exists a basis, it can be traded (basis swap).
Implied forward rate $$\tilde{L_T} = \frac{C_T \sum_{t=1}^T DF_t - \sum_{t=1}^{T-1} L_t DF_t}{DF_T} \equiv \frac{CD - LD}{DF_T}$$
Note that
$$\begin{aligned} PV(Fix) &=PV(Var)\\ C_T \sum_{t=1}^T DF_t &= \sum_{t=1}^T \tilde{L_t} DF_t = \sum_{t=1}^{T-1} \tilde{L_t} DF_t + \tilde{L_T} DF_T\\ \end{aligned}$$"""Calculate implied forward rates (multi curve approach)
"""
def multi_curve(swap_rates, eonia_rates):
# bootstrap discount factors from eonia rates
discount_factors = swap_bootstrap(eonia_rates)[0]
# setup
SDF = 0 # sum of discount factors
CD = 0 # discounted fixed leg until T
LD = 0 # discounted float leg until T-1
implied_forwards = []
for C_T, DF_T in zip(swap_rates, discount_factors):
SDF += DF_T # add up discount factors
CD = C_T * SDF # C_T is an annuity
L = (CD - LD) / DF_T # calculate implied forward rate (see equation above)
implied_forwards.append(round(L,4)) # round and save L
LD += L * DF_T # sum up discounted L (used in next iteration)
return implied_forwards
print(f"Implied forward rates: {multi_curve([0.04, 0.05, 0.06], [0.033, 0.045, 0.053])}")
Implied forward rates: [0.04, 0.0606, 0.082]
"""Price a forward swap (multi curve approach)
"""
implied_forwards = multi_curve([0.04, 0.05, 0.06], [0.033, 0.045, 0.053])
eonia_discount_factors = swap_bootstrap([0.033, 0.045, 0.053])[0]
PV_Var = sum(np.array(implied_forwards[1::]) * np.array(eonia_discount_factors[1::]))
s_13m = PV_Var / sum(eonia_discount_factors[1::])
print(f"1X2 swap rate (multi curve) = {s_13m*100:.2f}%")
1X2 swap rate (multi curve) = 7.09%
CCY basis swap:
In theory the spread of the product should be almost zero. However due to different liquidity levels and different credit levels in the currencies in reality we find substantial spreads.
Let us start with two definitions:
The FX swap is an agreement to exchange spot at $S$ as well as in the future at forward rate $F_T$. According to Covered Interest Rate Parity, this has arbitrage relationship:
$$ \begin{align*} (1 + r_B \times T_B) F_T &= S (1 + r_Q \times T_Q) \\ &...\\ F_T &= S \frac{(1 + r_Q \times T_Q)}{(1 + r_B \times T_B)} \end{align*} $$By convention, FX forwards (aka FX ourights) are often quotet in FX swap points, this is,
$$\text{FX swap points} = F_T - S$$The FX swap market includes besides the interest rate differential also liquidity and credit advantages of the two currencies. Since the financial crisis the interest differential (interest rate parity) is not a good indicator for the FX swap anymore.
A Forward Rate Agreement (or FRA) is an agreement between two parties to exchange a short term underlying interest rate in the future.
$$\text{Settlement amount} = \frac{(L(d_2)-r(d_3,d_4)) \frac{act}{360}}{\left(1+L(d_2) \frac{act}{360}\right)} \times \pm \text{Nominal}$$A buyer enters into the contract to protect himself from increasing interest rates, hence, $+ \rightarrow$ buyer, $- \rightarrow$ seller.
We will use a Bund future to highlight the peculiarities of future clearing, settlement and hedging. The Bund future is a future traded at EUREX on a hypothetical long term German government bond.
A clearing mechanism is a margining mechanism the involves one or more of:
Bund future specifications:
The ideal bond is not available, thus, we need a conversion factor $\text{CF}_{\text{Bond}}$
$$ \begin{align*} \text{CF}_{\text{Bond}} &= \frac{\text{PV}_{\text{Bond}}(@6\%)}{\text{PV}_{\text{Underlying}}(@6\%)}\\ &= \frac{\text{PV}_{\text{Bond}}(@6\%)}{100}. \end{align*} $$The conversion factor allows us to calculate the payment at deilvery
$$\text{Delivery payment} = \text{EDSP} \times \text{CF}_{\text{Bond}} \times 100,000 + \text{Accrued} $$“The final settlement price [EDSP] is established by Eurex on the final settlement day at 12:30 CET based on the volume-weighted average price of all trades during the final minute of trading provided that more than ten trades occurred during this minute; otherwise the volume-weighted average price of the last ten trades of the day, provided that these are not older than 30 minutes. If such a price cannot be determined, or does not reasonably reflect the prevailing market conditions, Eurex will establish the final settlement price.” (Eurex Exchange 2019)
Now that we have the conversion factor ($\text{CF}_{\text{Bond}}$), we have to choose which bond to deliver at maturity. This is the cheapest-to-deliver, which can be identified by its implied repo rate
$$ \text{Implied repo} = \frac{\text{Delivery payment} - \text{Dirty price}}{\text{Dirty price} \times \frac{\text{Days to settlement}}{360}} $$A Repo trade is the agreement to sell (Spot) and buy back ($t+n$) an asset. Spot short and at the same time future long can be seen as a repo trade. Another interpretation is collateralized lending. The bond with the highest Implied Repo rate is the cheapest to deliver (CTD).
The forward price for the bond cheapest to deliver (CTD) is determined by the costs of acquiring the bond now and holding it until it shall be delivered to fulfil the obligation from the future. The cost of carry model applies, hence,
$$\text{Forward price} = \text{Spot} + \text{Financing} - \text{Accrued}$$Considering the underlying of the future,
$$\text{Theoretical future} = \frac{\text{Forward price}}{\text{CF}_{\text{Bond}}}$$Basis Point Value (BPV = PVO1) Hedge is a Nominal hedge, adjusted for underlying Basis Risk.
$$\text{Basis point hedge ratio} = \text{Simple hedge ratio} \times \text{CF}_{\text{Bond}} \times \frac{\text{BPV}_\text{Spot}}{\text{BPV}_\text{CTD}}$$PVO1 or BPV is the MTM change caused by a $0.0001$ = 1bp shift in yield.
Notation:
| Lower boundaries | Upper boundaries | |
|---|---|---|
| European Call | $c \ge S-K e^{-rT}$ | $c \le S$ |
| American Call | $C \ge S-Ke^{-rT}$ | $C \le S$ |
| European Put | $p \ge K e^{-rT} - S$ | $p \le K e^{-rT}$ |
| American Put | $P \ge K - S$ | $P \le K$ |
Why it makes no sense to execute an american call option prematurely:
Consider random variable $X \sim \mathcal{N}(\mu,\sigma^2)$ as an example for an absolutely continuous univariate distribution.
Probability density function (PDF): $$\phi(x) = \frac{d\Phi(x)}{dx} = \frac{1}{\sigma \sqrt{2\pi}} e^{-\frac{1}{2} \left( \frac{x-\mu}{\sigma} \right)^2 }$$
Cumulative distribution function (CDF): $$\Phi(x) = P(X \le x) = \int_{-\infty}^x \phi(u)\; du$$
Quantile function, i.e. inverse of the CDF: $$\Phi^{-1}(p) = \inf \{x \in \mathbb{R}: p \le \Phi(x)\},\;\; p \in (0,1)$$
"""PDF vs CDF
"""
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import norm
x = np.linspace(-3, 3, 100)
p = np.linspace(0, 1)
plt.figure("pdf")
plt.title("Probability Density Function (PDF)")
plt.plot(x, norm.pdf(x, 0, 1))
plt.plot(x, norm.pdf(x, 1, 1))
plt.plot(x, norm.pdf(x, 0, 1.5))
plt.axvline(x=-2.5, color='red', linewidth=.75)
plt.axvline(x=0, color='grey', linewidth=.75)
plt.legend(["N(0, 1)", "N(1, 1)", "N(0, 1.5)"], loc='upper left')
plt.figure("cdf")
plt.title("Cumulative Distribution Function (CDF)")
plt.plot(x, norm.cdf(x, 0, 1))
plt.plot(x, norm.cdf(x, 1, 1))
plt.plot(x, norm.cdf(x, 0, 1.5))
plt.axvline(x=-2.5, color='red', linewidth=.75)
plt.axvline(x=0, color='grey', linewidth=.75)
plt.legend(["N(0, 1)", "N(1, 1)", "N(0, 1.5)"], loc='upper left')
plt.figure("qf")
plt.title("Quantile function (inverse CDF)")
plt.plot(p, norm.ppf(p, 0, 1))
plt.plot(p, norm.ppf(p, 1, 1))
plt.plot(p, norm.ppf(p, 0, 1.5))
plt.axvline(x=0.5, color='grey', linewidth=.75)
plt.legend(["N(0, 1)", "N(1, 1)", "N(0, 1.5)"], loc='upper left')
plt.show()
In probability theory and statistics, the term Markov property refers to the memoryless property of a stochastic process, i.e. only the present price is relevant for predicting the future. Markov property is consistent with the weak form of market efficiency!
A stochastic process $X$, with drift $\mu$ and dispersion $\sigma$, follows a generalized Wiener process if it satisfies the following stochastic differential equation (SDE)
$$dX_t = \mu \times dt + \sigma \times dW_t$$Note that by definition, the increments of a wiener process $W_t = W_t - W_0 \sim \mathcal{N}(0,t)$, i.e. are normally distributed, centered at zero.
For an arbitrary initial value $X_0$
$$X_t = X_0 + \mu \times t + \sigma \times z \times \sqrt{t}$$, where $z \sim \mathcal{N}(0,1)$.
"""Calculate the probability of illiquidity
"""
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import norm
month = range(44)
def prob_illiquid(m, const=2, drift=.1, variance=.16):
expected_cash = const + drift * m
sigma = (variance * m)**0.5
p_illiquid = norm.cdf( (0 - expected_cash) / sigma ,0,1) # normalization: (x - mu) / sigma
return p_illiquid
print(f"Probability of illiquidity in 6 month: {prob_illiquid(6)*100:.4f}%")
print(f"Probability of illiquidity in 12 month: {prob_illiquid(12)*100:.4f}%")
print(f"Probability of illiquidity in 36 month: {prob_illiquid(36)*100:.4f}%")
X = range(1,12*5, 1)
fX = list(map(lambda m: prob_illiquid(m), X))
plt.figure('Illiquidity', figsize=(9,4))
plt.title("Illiquidity")
plt.plot(X, fX)
plt.xlabel("Month (m)")
plt.ylabel("Probability")
plt.show()
Probability of illiquidity in 6 month: 0.3982% Probability of illiquidity in 12 month: 1.0461% Probability of illiquidity in 36 month: 0.9815%
A stochastic process $S$ follows a Geometric Brownian Motion (GBM) if it satisfies the following stochastic differential equation (SDE)
$$dS_t = \mu S_t dt + \sigma S_t dW_t$$where $W_t = W_t - W_0 \sim \mathcal{N}(0,t)$, i.e. a Wiener process, hence, normally distributed increments.
For an arbitrary initial value $S_0$, the above SDE has the following analytical solution (under Ito's interpretation)
$$S_t = S_0 \exp\left( \left(\mu - \frac{\sigma^2}{2}\right) t + \sigma W_t \right)$$, hence,
$$\ln \frac{S_t}{S_0} = \left(\mu - \frac{\sigma^2}{2}\right) t + \sigma W_t$$Note that $W_t = z \sqrt{t}$, where $z \sim \mathcal{N}(0,1)$. Therefore, $\ln \frac{S_t}{S_0}$ (log-return) is normally and $S_t$ is log-normally distributed.
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
S_0 = 50
mu = 0.16
sigma = 0.3
# Expected price and price volatility
S_1 = 50 * (1 + mu / 365)
print(f"E[S] tomorrow = {S_1:.4f}")
price_vola = 50 * sigma * (1 / 365)**0.5
print(f"Price vola = {price_vola:.4f}")
# Confidence interval
q975 = norm.ppf(0.975, 0, 1)
print(f"With 95% confidence: {S_1 - price_vola * q975 :.2f} < S_1 < {S_1 + price_vola * q975 :.2f}")
# Value at risk
q95 = norm.ppf(0.95, 0, 1)
print(f"1 day 95% VaR = {S_0 - (S_1 + price_vola * q95) :.2f}")
def gbm(mu, sigma, N=15):
"""Cumulative returns from a GBM
"""
# Random numbers N(0,1)
Z = np.random.normal(0, 1, N)
# returns for our geometric brownian motion
R = Z * sigma * np.sqrt(1/250) + (mu - 0.5*sigma**2) * 1/250
# stock process (insert 100 for t=0)
S = np.insert(np.exp(np.cumsum(R)), 0, 1)
return S
plt.figure('MC', figsize=(9,4))
plt.title("Price path")
for n in range(100):
plt.plot(S_0*gbm(mu, sigma), linewidth=2)
plt.xlabel("Trading days")
plt.ylabel("Price")
plt.show()
E[S] tomorrow = 50.0219 Price vola = 0.7851 With 95% confidence: 48.48 < S_1 < 51.56 1 day 95% VaR = -1.31
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets
import matplotlib.pyplot as plt
%matplotlib inline
#%% Create stock process:
# Random numbers N(0,1)
np.random.seed(888)
Z = np.random.normal(0, 1, 250)
def plot_gbm(mu, sigma):
# returns for our geometric brownian motion
R = Z * sigma * np.sqrt(1/250) + (mu-0.5*sigma**2) * 1/250
# stock process (insert 100 for t=0)
S = np.insert(100*np.exp(np.cumsum(R)), 0, 100)
# Plotting
plt.figure('GBM process', figsize=(9,4))
plt.title(f"GBM ($\mu$ = {mu}, $\sigma$ = {sigma})")
plt.axhline(y = 100, color ="black", linestyle ="--", linewidth=1)
plt.plot(S)
plt.xlabel("Time")
plt.ylabel("Price")
plt.ylim(30, 170)
plt.show()
interact(plot_gbm,
mu = widgets.FloatSlider(value=0.1, min=-0.5, max=0.5, step=0.05),
sigma = widgets.FloatSlider(value=0.25, min=0, max=0.9, step=0.05),
)
<function __main__.plot_gbm(mu, sigma)>
"""Monte Caro simulation vs Black Scholes
"""
import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt
# Option Parameters
S_0 = 220 # price of the underlying
K = 220 # strike price
sigma = .98 # volatility
rate = np.log(1.21) # converting the ISDA rate to continous
T = 1 # time to maturity
steps = 2 # steps
# We talk about this later, here its just a reference point
def black_scholes(cpflag,S,K,T,r,sigma):
# cpflag in {1 for call, -1 for put}
d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
d2 = (np.log(S / K) + (r - 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
price = cpflag * (S * norm.cdf(cpflag*d1, 0.0, 1.0) - K * np.exp(-r * T) *
norm.cdf(cpflag*d2, 0.0, 1.0))
return price
# Monte carlo pricing
def monte_carlo(cpflag,S,K,T,r,sigma,Nsim):
Z = np.random.normal(0, 1, Nsim)
S_T = S * np.exp(Z * sigma * np.sqrt(T) + (r-0.5*sigma**2) * T)
CF_T = np.maximum(cpflag*(S_T - K),0)
price = np.mean(CF_T) * np.exp(-r*T)
return price
# Increase Nsim from 1000 to 100000 by 1000
X = range(1000,100000, 1000)
# Compute prices for all Nsim in X
mc_prices = list(map(lambda n: monte_carlo(1,220,220,1,np.log(1.21),0.98,n), X))
plt.figure('Monte Caro simulation vs Black Scholes', figsize=(9,4))
plt.plot(X,mc_prices)
plt.axhline(y=black_scholes(1,220,220,1,np.log(1.21),0.98), color='orange')
plt.xlabel("Number of simulations")
plt.ylabel("Price of the option")
plt.legend(["Monte Carlo", "Black-Scholes"])
plt.show()
black_scholes(1,220,220,1,np.log(1.21),0.98)
95.91455141466928
Consider a replicating portfolio with bond (cash) $B$ and $\Delta$ share in the underlying asset $S$
$$ \begin{align} V_0 &= B + \Delta S_0\\ V_u &= B e^{rt} + \Delta S_u = (S_u - K)^+\\ V_d &= B e^{rt} + \Delta S_d = (S_d - K)^+ \end{align} $$Note that
$$\Delta = \frac{(S_u - K)^+ - (S_d - K)^+}{S_u - S_d}$$"""Replicate the option payout in a one step lattice
"""
import numpy as np
def delta(flag, S_u, S_d, K): # flag = 1 (Call); -1 (Put)
delta = ( np.maximum(flag*(S_u - K), 0) - np.maximum(flag*(S_d - K), 0) ) / (S_u - S_d)
return delta
#print(delta(1, 60, 40, 45))
def bond(flag, S_u, S_d, K, r, T):
beta = (np.maximum(flag*(S_u - K), 0) - delta(flag, S_u, S_d, K) * S_u) * np.exp(-r*T)
return beta
#print(beta(1, 60, 40, 45, 0.05, 1))
print(f"V_u = {bond(1, 60, 40, 45, 0.05, 1)*np.exp(.05)+delta(1, 60, 40, 45)*60:.0f}")
print(f"V_d = {bond(1, 60, 40, 45, 0.05, 1)*np.exp(.05)+delta(1, 60, 40, 45)*40:.0f}")
print(f"V_0 = {bond(1, 60, 40, 45, 0.05, 1)+delta(1, 60, 40, 45)*50:.4f}")
V_u = 15 V_d = -0 V_0 = 8.9631
The underlying of an option with maturity $T$ is modelled by an $n$ step binomial tree, with time between two steps of the tree $\Delta t = T / n$. With up move $u$ and down move $d$
$$u = \exp\left(\sigma \sqrt{\Delta t}\right)$$$$d = u^{-1}$$The following must hold for the up move probability $p$, in order to ensure martingale property,
$$S_0 = \left( S_0 \times u \times p + S_0 \times d \times (1-p) \right) e^{-r\Delta t} \implies p = \frac{e^{r\Delta t} - d}{u - d}$$Remember the binomial probability mass function (PMF)
$$f(k,n,p) = {n \choose k} p^k (1-p)^{n-k}$$here, $n$ is the number of steps in the tree and $k$ the number of up moves to reach a state $S_T$.
The CRR model is solved recuresively by discounting expected values. For American options, the buyer has the right to exercise at any point in time. Due to its discrete modelling of time, the CRR model allows execution at any node in the lattice. However, the buyer would only execute if the value that can be realized through execution $IV_t$ (intrinsic value) is above the discounted expected value of the next node. Therefore, the value that is considered for each node is
$$V_t = \max\left(E[V_{t+\Delta t}]e^{-r\Delta t}, IV_t \right).$$Solve recursively by stepwise discounting $V_t$ backwards in the tree. In this example we consider a European call with $S_0 = 220$, $K=165$, $\sigma = 98\%$, $T=1$, $r = 21\%$ ISDA.
"""Replicating options, Cox Ross Rubinstein style
"""
import numpy as np
import matplotlib.pyplot as plt
# INPUTS
S_0 = 220 # price of the underlying
K = 220 # strike price
sigma = .98 # volatility
rate = np.log(1.21) # converting the ISDA rate
T = 1 # time to maturity
steps = 2 # steps
# see Slide 10
dt = T / steps
u = np.exp(sigma * np.sqrt(dt))
d = 1 / u
binomial_tree = []
for n in range(steps+1):
binomial_tree.append([])
for k in range(n+1):
S_nk = S_0 * u**(k) * d**(n-k)
binomial_tree[-1].append(S_nk)
plt.plot(n*dt, S_nk, 'ro') # this will take forever for large # steps
plt.annotate(f"{S_nk:.0f} ", (n*dt, S_nk), ha='right', va='center')
print(binomial_tree)
plt.show()
[[220.0], [110.02008068234205, 439.9196919310029], [55.02008251522298, 220.0, 879.6787970394022]]
"""Full CRR model
"""
import matplotlib.pyplot as plt
import numpy as np
def Binomial(n, S, K, r, v, t, PutCall=1, American=False):
At = t/n
u = np.exp(v*np.sqrt(At))
d = 1./u
p = (np.exp(r*At)-d) / (u-d)
# Binomial price tree
stockvalue = np.zeros((n+1,n+1))
stockvalue[0,0] = S
for i in range(1,n+1):
stockvalue[i,0] = stockvalue[i-1,0]*u
for j in range(1,i+1):
stockvalue[i,j] = stockvalue[i-1,j-1]*d
# option value at final node
optionvalue = np.zeros((n+1,n+1))
for j in range(n+1):
optionvalue[n,j] = max(0, PutCall*(stockvalue[n,j]-K))
# recursive calculations
for i in range(n-1,-1,-1):
for j in range(i+1):
if American:
optionvalue[i,j] = max(0, PutCall*(stockvalue[i,j]-K), np.exp(-r*At)*(p*optionvalue[i+1,j]+(1-p)*optionvalue[i+1,j+1]))
else:
optionvalue[i,j] = np.exp(-r*At)*(p*optionvalue[i+1,j]+(1-p)*optionvalue[i+1,j+1])
return optionvalue[0,0]
# Inputs
S_0 = 220 # price of the underlying
K = 220 # strike price (165 vs 220)
sigma = .98 # volatility
rate = np.log(1 + 0.21) # converting the ISDA rate
T = 1 # time to maturity
steps = 2 # steps
#Graphs and results for the Option prices
print("European Call Price: %s" %(Binomial(steps, S_0, K, rate, sigma, T, PutCall=1)))
print("American Call Price: %s" %(Binomial(steps, S_0, K, rate, sigma, T, PutCall=1, American=True)))
print("European Put Price: %s" %(Binomial(steps, S_0, K, rate, sigma, T, PutCall=-1)))
print("American Put Price: %s" %(Binomial(steps, S_0, K, rate, sigma, T, PutCall=-1, American=True)))
European Call Price: 87.25679654189224 American Call Price: 87.25679654189224 European Put Price: 49.07497836007403 American Put Price: 59.982962882003555
from scipy.stats import binom
# Recall our CRR parameters
sigma = .98 # volatility
rate = np.log(1 + 0.21) # converting the ISDA rate
T = 1 # time to maturity
steps = 2 # steps
dt = T / steps
u = np.exp(sigma*np.sqrt(dt))
d = 1 / u
p = (np.exp(rate*dt)-d) / (u-d)
CoN_payout = []
AoN_payout = []
for n, S in enumerate(binomial_tree[-1]): # This is the binomial tree that we generated above
if S > 165: # strike
AoN_payout.append( binom.pmf(n, steps, p) * S * np.exp(-rate*T) )
CoN_payout.append( binom.pmf(n, steps, p) * 165 * np.exp(-rate*T) )
AoN = sum(AoN_payout)
CoN = sum(CoN_payout)
print(f"Price asset or nothing = {AoN:.4f}")
print(f"Price cash or nothing = {CoN:.4f}")
print(f"AoN - CoN = {AoN - CoN:.4f}")
print(f"CRR call price (above) = {Binomial(steps, S_0, 165, rate, sigma, T, PutCall=1):.4f}")
Price asset or nothing = 203.6337 Price cash or nothing = 87.2827 AoN - CoN = 116.3510 CRR call price (above) = 116.3510
with
$$d_\pm = \frac{\ln\left(\frac{S_0}{K}\right) + rT \pm \sigma^2 \frac{T}{2}}{\sigma \sqrt{T}}$$with
$$d_\pm = \frac{\ln\left(\frac{S_0}{K}\right) + (r-q)T \pm \sigma^2 \frac{T}{2}}{\sigma \sqrt{T}}$$with
$$d_\pm = \frac{\ln\left(\frac{F_T}{K}\right) \pm \sigma^2 \frac{T}{2}}{\sigma \sqrt{T}}$$"""Black Scholes sensitivities
"""
import numpy as np
from scipy.stats import norm
#%% Black-Scholes formula
def black_scholes(cpflag,S,K,T,r,sigma):
# cpflag in {1 for call, -1 for put}
d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
d2 = (np.log(S / K) + (r - 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
price = cpflag * (S * norm.cdf(cpflag*d1, 0.0, 1.0) - K * np.exp(-r * T) *
norm.cdf(cpflag*d2, 0.0, 1.0))
return price
print("European Call Price: %s" %(black_scholes(1,220,220,1,np.log(1.21),.98)))
print("European put Price: %s" %(black_scholes(-1,220,220,1,np.log(1.21),.98)))
European Call Price: 95.91455141466928 European put Price: 57.73273323285109
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets
import matplotlib.pyplot as plt
%matplotlib inline
def plot_block_scholes(pcflag, sigma, r, T):
# Initial parameters
X = np.linspace(1,100,100)
plt.figure('BSM Price', figsize=(9,4))
plt.plot(X,np.maximum(pcflag*(X-50),0), '--')
plt.plot(X,black_scholes(pcflag,X,50,1,0.05,0.25))
plt.plot(X,black_scholes(pcflag,X,50,T,r,sigma))
plt.legend(["Intrinsic value", "Initial parameters", "Modified parameters"])
plt.show()
interact(plot_block_scholes,
pcflag=widgets.FloatSlider(value=1, min=-1, max=1, step=2),
sigma=widgets.FloatSlider(value=0.25, min=0.1, max=0.5, step=0.05),
r=widgets.FloatSlider(value=0.05, min=0.01, max=0.1, step=0.01),
T=widgets.FloatSlider(value=1, min=0.5, max=5, step=0.5)
)
<function __main__.plot_block_scholes(pcflag, sigma, r, T)>
"""CRR vs Black Scholes
"""
print("European Call Price: %s" %(black_scholes(1,220,220,1,np.log(1.21),.98)))
print("European put Price: %s" %(black_scholes(-1,220,220,1,np.log(1.21),.98)))
X = range(2,100, 1)
crr_prices = list(map(lambda n: Binomial(n, 220, 220, np.log(1.21), .98, 1, PutCall=1), X))
plt.figure('CRR vs Black Scholes', figsize=(9,4))
plt.plot(X,crr_prices)
plt.axhline(y=black_scholes(1,220,220,1,np.log(1.21),.98), color='orange')
plt.xlabel("Number of steps in the tree")
plt.ylabel("Price of the option")
plt.legend(["Cox-Ross-Rubinstein", "Black-Scholes"])
plt.show()
European Call Price: 95.91455141466928 European put Price: 57.73273323285109
The European call price is simply the discounted expected value of its payout $E^\mathbb{Q}[(S_T - K)^+] e^{-rT}$, with
$$\begin{align} E^\mathbb{Q}[(S_T - K)^+] &= \int_K^\infty (S_T - K) f(S_T) dS_T\\ &=\int_K^\infty S_T f(S_T) dS_T - K \int_K^\infty f(S_T) dS_T \end{align}$$The terminal stock price at time $T$, under the risk-neutral measure, follows a log-normal distribution with $\mu = \ln S_0 + (r_f-\sigma^2 / 2)T$, variance $s^2 = \sigma^2 T$ and PDF
$$f(x) = \frac{1}{x s \sqrt{2 \pi}} \exp\left( -\frac{(\ln x - \mu)^2}{2 s^2} \right) $$Consider a contract where you receive the underlying $S_T$ conditional on $S_T > K$, note that this is different from a call option where the cash flow is $S_T - K$ under the condition that $S_T > K$. The price of this (expensive) option is $E^\mathbb{Q}[S_T 1_{\{S_T > K\}}] e^{-rT}$, where
$$E^\mathbb{Q}[S_T 1_{\{S_T > K\}}] = \int_K^\infty S_T f(S_T) dS_T$$Now consider an option that pays some amount $K$ (cash) conditional on $S_T > K$. The price of this (again expensive) option is $K \times P^\mathbb{Q}[S_T > K] e^{-rT}$, where
$$P^\mathbb{Q}[S_T > K] = \int_K^\infty f(S_T) dS_T$$Therefore, a European call option is a portfolio out of a long asset-or-noting option and a short cash-or-nothing option.
"""Straightforward integration vs Black Scholes
"""
import numpy as np
from scipy.stats import norm
from scipy.integrate import quad
import matplotlib.pyplot as plt
# See Chapter 2.4
def black_scholes(cpflag,S,K,T,r,sigma):
# cpflag in {1 for call, -1 for put}
d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
d2 = (np.log(S / K) + (r - 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
price = cpflag * (S * norm.cdf(cpflag*d1, 0.0, 1.0) - K * np.exp(-r * T) *
norm.cdf(cpflag*d2, 0.0, 1.0))
return price
# probability density function P(S_T > K)
def PdF(S, S_0, rf, sd, t):
A = 1 / ( np.sqrt(2*np.pi*t)*sd*S)
B = np.exp(-(np.log(S)-(np.log(S_0)+(rf-0.5*sd**2)*t))**2/(2*sd**2*t))
return A * B
# expected value S_T
def EdF(S, S_0, rf, sd, t):
A = 1 / ( np.sqrt(2*np.pi*t)*sd*S)
B = np.exp(-(np.log(S)-(np.log(S_0)+(rf-0.5*sd**2)*t))**2/(2*sd**2*t))
return S * A * B
# expected value of (S_T - K)
def E2dF(S, S_0, rf, sd, t, K):
A = 1 / ( np.sqrt(2*np.pi*t)*sd*S)
B = np.exp(-(np.log(S)-(np.log(S_0)+(rf-0.5*sd**2)*t))**2/(2*sd**2*t))
return (S - K) * A * B
# Inputs
S_0 = 220 # price of the underlying
K = 220 # strike price
sigma = .98 # volatility
rate = np.log(1.21) # converting the ISDA rate
T = 1 # time to maturity
steps = 2 # steps
# asset or nothin
AoN = quad(lambda x: EdF(x, S_0, rate, sigma, T), K, +np.inf)[0] * np.exp(-rate*T)
# cash or nothing with cash = K
CoN = K * quad(lambda x: PdF(x, S_0, rate, sigma, T), K, +np.inf)[0] * np.exp(-rate*T)
# straight forward integration
Call = quad(lambda x: E2dF(x, S_0, rate, sigma, T, K), K, +np.inf)[0] * np.exp(-rate*T)
print(f"Black-Scholes price : {black_scholes(1,S_0,K,T,rate,sigma):.8f}")
print(f"AoN - CoN : {AoN-CoN:.8f}")
print(f"Straightforward integration: {Call:.8f}")
x = np.linspace(1, 800, 100)
fig, ax1 = plt.subplots(figsize=(9,4))
ax1.plot(x, PdF(x, S_0, rate, sigma, T))
ax1.set_xlabel('S_T')
ax2 = ax1.twinx()
ax2.plot(x, np.maximum(x-K,0), '--', color="black")
ax2.set_ylabel('max(S_T - K, 0)')
plt.show()
Black-Scholes price : 95.91455141 AoN - CoN : 95.91455141 Straightforward integration: 95.91455141
Greeks are partial derivatives that describe the sensitivity of the price of an option to a change in the respective underlying parameters.
Greeks are not constant, see e.g. delta (red) and gamma (blue)
All parameters that determine the option price in a BSM model are observable, except $\sigma$. However, option prices are not calculated, they are the result of supply and demand on a liquid market. Thus, given an option price, we can find the $\sigma$ that ceteris paribus equates the BSM price with the market price.
"""Calculate implied volatility
"""
# Inputs
S_0 = 220 # price of the underlying
K = 220 # strike price
rate = np.log(1.21) # converting the ISDA rate
T = 1 # time to maturity
steps = 2 # steps
P_MKT = 85 # market price of the option
def iv_solve(IV):
P_BSM = black_scholes(1,S_0,K,T,rate,IV)
return (P_BSM - P_MKT)**2
IV = minimize(iv_solve, x0=0.25, bounds = ((0, None),)).x[0]
print(f'Implied vol: {IV:.4f}')
print(f'BSM price : {black_scholes(1,S_0,K,T,rate,IV):.4f}')
Implied vol: 0.8251 BSM price : 85.0000
Consider a matrix of option prices for different strikes and maturities. According to the theory, the resulting volatility surface should be flat, however, we usually observe sth like
Reasons for a non-flat volatility surface
A CDS contract works like an insurence that is linked to a credit event. This is, the protection buyer pays a premium $S$ (CDS spread) to the protection seller until the contract matures in $T$. In case of a default during the term of the contract, i.e. default in $\tau < T$, the protection buyer receives a settlement payment from the seller. This settlement payment is often, but not necessarily, designed to cover the loss given the default of a credit risky position, i.e. $(1 - \mathrm{Recovery})$.
Credit events (ISDA credit derivative definition 5/99 6/03 4/09 2/14):
The payout can be replicarted with one credit risky and another risk free floater, which both trade at par. Shorting the risky floater yields 100% nominal, which can be invested into the risk free floater. The latter pays per definition the risk free rate $r_f$, wheres payments into the risky floater include an additional (credit) spread $r_f + S_{FRN = Par}$. In case of a credit event, the risk free floater is sold at 100%, whereas the risky floater is bought back at recovery. As this replicates the payout of the CDS $(1 - \mathrm{Recovery})$, the replicating credit spread must be equal to the CDS spread $S$. However, this is only an approximation, as it ignores that recovery is only paid on the nominal amount, transaction costs, and liquidity differences (availability) of the floaters.
$$S_{FRN = Par} \approx S$$If a floater is not trading at par, the replication does not work. However, the PV difference must be due to the discounted difference in spreads. Thus, given the spread $S_{FRN \neq Par}$ and price $P_{FRN\neq Par}$ of a non par floater, we can solve for the spread of the corresponding par floater $S_{FRN = Par}$
$$ \begin{align} P_{FRN\neq Par} - 100\% &= PV(S_{FRN\neq Par} - S_{FRN = Par})\\ &= (S_{FRN\neq Par} - S_{FRN = Par}) \sum_{t=1}^T DF_t\\ &...\\ S_{FRN = Par} &= S_{FRN\neq Par} - \frac{P_{FRN \neq Par} - 100\%}{\sum_{t=1}^T DF_t} \end{align} $$In reality, costs arise from shorting the FRN, i.e. the bond has to be borrowd for a fee $S_{Repo}$. Thus
$$S_{FRN = Par} + S_{Repo} = S$$No floater available? No problem, use an asset swap to create a synthetic floater, which at par is worth
$$ \begin{align} 100\% &= P_\mathrm{Bond} - PV(C_\mathrm{Bond}) + PV(\mathrm{Libor} + S)\\\\ &= P_\mathrm{Bond} - C_\mathrm{Bond} \sum_{t=1}^T DF_t + 1 - DF_T + S \sum_{t=1}^T DF_t\\\\ S &= \frac{DF_T + C_\mathrm{Bond} \sum_{t=1}^T DF_t - P_\mathrm{Bond}}{\sum_{t=1}^T DF_t} \end{align} $$If the bond is trading at par, $S$ is simply
$$S = C_\mathrm{Bond} - \frac{1-DF_T}{\sum_{t=1}^T DF_t}\\\\$$The asset swap spread is only an approximation for the CDS spread, as the CDS cash flow is not perfectly replicated (close out value of the swap) and there might be liquidity issues.
As always, the price of a financial product is the sum of its discounted (expected) cash flows. Accounting for credit risk simply means accounting for both, (risk free) interest rate and the probability of default.
Let $h_t$, the hazard rate (default intensity), be the probability to default in an infinitesimally short time interval $dt$.
$$h_t = \ln\left(\frac{1-p_{t-1}}{1-p_t} \right) \implies e^{h_t} = \frac{1-p_{t-1}}{1-p_t} \Longleftrightarrow p_t = 1 - \frac{1- p_{t-1}}{e^{h_t}}$$Then $a_t(h_t)$ is the present value of a paymment of $1$ at time $t$ if $t < \tau$, hence, no default. Think of this as a discount factor ($DF_t$) that also accounts for the likelyhood of the asset being alive.
$$a_t(h_t) = e^{-(h_t+r_t)t}$$On the other hand, $b_t(h_t)$ is the present value of a payment of 1 at time $t$ if $t-1 < \tau < t$, hence, default between $t-1$ and $t$.
$$b_t(h_t) = e^{-r_t t} \left( e^{-(h_t)t_{-1}} - e^{-(h_t)t} \right)$$Summing up the risky discount factors $a_t(h_t)$ and $b_t(h_t)$
$$ \begin{align} A_T &= \sum_{t=1}^T a_t\\ B_T &= \sum_{t=1}^T b_t \end{align} $$Allows us to price a CDS with spread $S$ as
$$ \begin{align} 0 &= PV(1-\mathrm{Recovery}) - PV(S)\\ &= B_T(1-\mathrm{Recovery}) - A_T S\\\\ S &= \frac{B_T(1-\mathrm{Recovery})}{A_T} \end{align} $$Hint: compare this to the single curve swap pricing above.
In practise, in preparation for central counterparty clearing, spread payments into a CDS are fixed at, i.e. $S^* \in \{0.002,0.01,0.05,0.1\}$, see ISDA agreements. Also the recovery rate is fixed, e.g. to 40% for investment grade. Therefore, contracts are unfair when entered and require an upfont payment
$$\mathrm{Upfront} = B_T(1-\mathrm{Recovery}) - A_T S^*$$"""Fair spread and upfront payment
"""
import numpy as np
def CDS_AB(h, r, T):
# assuming a constant hazard rate h
A_T = 0
B_T = 0
for t in range(1, T+1):
r_t = r[t-1] if isinstance(r, (list, np.ndarray)) else r
A_T += np.exp(-(h + r_t)*t)
B_T += np.exp(-r_t*t) * ( np.exp(-h*(t-1)) - np.exp(-h*t) )
return (A_T, B_T)
# Example:
h = 0.0209224
R = 0.3
# Fair spread:
A_T, B_T = CDS_AB(h, np.log([1.0365, 1.0386]), 2)
S = B_T*(1 - R) / A_T
print(f"Fair CDS spread: {S:.4f}")
# Upfront payment (spread payment of 100 basis points = 0.01)
upfront = (B_T*(1 - R) - A_T * 0.01) *10**7
upfront_alt = (S - 0.01) * A_T *10**7
print(f"Upfront payment on 10mln nominal: {upfront:.4f}")
print(f"Upfront payment on 10mln nominal: {upfront_alt:.4f}")
Fair CDS spread: 0.0148 Upfront payment on 10mln nominal: 88025.0332 Upfront payment on 10mln nominal: 88025.0332
Recall the CDS replication via floating rate notes, solve for $h$ in $A_T$ and $B_T$, using that $S_{FRN = Par} \approx S$
$$ \begin{align} P_{FRN} - 100\% &= (S_{FRN} - S) A_T\\ &= A_T S_{FRN} - B_T(1-\mathrm{Recovery})\\ \end{align} $$Note that
$$\frac{S}{1-\mathrm{Recovery}} = \frac{B_T}{A_T} \approx h$$This approximation is often called "credit triangle".
"""Calibrate constant hazard rate to floater
"""
from scipy.optimize import minimize
r = 0.0917
T = 2
P_FRN = 1
S_FRN = 0.004
R = 0.2
def h_solve(h):
A_T, B_T = CDS_AB(h, r, T)
return (A_T*S_FRN - B_T*(1-R) + 1 - P_FRN)**2*100
# Find h that solves PV(alive) = PV(dead)
h = minimize(h_solve, x0=0.0001, bounds = ((0, None),)).x[0]
A_T, B_T = CDS_AB(h, r, T)
S = B_T*(1 - R) / A_T
print(f'Calibrated hazard rate: {h:.6f}')
print(f'Credit triangle : {S_FRN/(1-R):.6f}') # this works because P_FRN = 100% !
print(f'Calibrated CDS spread : {S:.6f}')
print(f'With A_T={A_T:.4f} and B_T={B_T:.6f}')
Calibrated hazard rate: 0.004988 Credit triangle : 0.005000 Calibrated CDS spread : 0.004000 With A_T=1.7320 and B_T=0.008660
# This will convert the notebook to a html file
import os
# Convert to html slides
os.system('jupyter nbconvert derivatives.ipynb --to html')
0